Clustering and Classification of Satellite Imagery

Abstract

-RAQUEL # Unsupervised Clustering -RAQUEL any intro stuff about this topic, feell free to use references This section will use unsupervised clustering to split up given satellite images into several groups, with the goal of identifying meaningfull clusters that represent different geographic features such as water, farm land, forest or buildings.

## Loading required package: randomForest
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'randomForest'
rasterJ<-brick("C:\\Users\\Owner\\Documents\\Portugal\\Sat_Image_Classification\\data\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
#rasterJ<-brick("C:\\Users\\D059331\\Desktop\\DM GIC\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
rasterJ
## class       : RasterBrick 
## dimensions  : 4238, 4239, 17964882, 4  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -105092.3, -94918.74, -50088.84, -39917.64  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\Owner\Documents\Portugal\Sat_Image_Classification\data\data\img\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF 
## names       : J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4 
## min values  :                                                  0,                                                  0,                                                  0,                                                  0 
## max values  :                                                255,                                                255,                                                255,                                                255
plotRGB(rasterJ, 3,2,1)

ext <- extent(-104294.4, -102964.5, -43623.48, -42742.44 )
sector <- crop(rasterJ, ext)
plotRGB(sector, 3, 2, 1)

sector[[5]] <- normdiff(sector)
plot(sector)

summary(sector)
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1
## Min.                                                    23
## 1st Qu.                                                 31
## Median                                                  34
## 3rd Qu.                                                 40
## Max.                                                   176
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2
## Min.                                                    26
## 1st Qu.                                                 40
## Median                                                  50
## 3rd Qu.                                                 65
## Max.                                                   255
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3
## Min.                                                     8
## 1st Qu.                                                 17
## Median                                                  34
## 3rd Qu.                                                 51
## Max.                                                   194
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4      values
## Min.                                                     4 -0.64885496
## 1st Qu.                                                 10 -0.20000000
## Median                                                  53  0.05645161
## 3rd Qu.                                                 70  0.15384615
## Max.                                                   191  0.76870748
## NA's                                                     0  0.00000000
means<-mov.fun(sector,3,mean, normalize=FALSE, verbose=TRUE)
## >> Band  1  of  5  done.
## >> Band  2  of  5  done.
## >> Band  3  of  5  done.
## >> Band  4  of  5  done.
## >> Band  5  of  5  done.
sds<-mov.fun(sector,3,sd, normalize=FALSE, verbose=TRUE)
## >> Band  1  of  5  done.
## >> Band  2  of  5  done.
## >> Band  3  of  5  done.
## >> Band  4  of  5  done.
## >> Band  5  of  5  done.
stk<-stack(means,sds)
sectorMovings<-brick(stk)
summary(sectorMovings)
##         layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1       layer.1
## Min.       24.00000    27.88889    9.111111    4.555556   -0.51877395
## 1st Qu.    30.88889    40.44444   16.666667   10.111111   -0.21001302
## Median     34.55556    51.00000   35.555556   53.555556    0.05963125
## 3rd Qu.    40.22222    64.88889   50.888889   70.222222    0.16173571
## Max.      172.11111   255.00000  189.222222  171.111111    0.73729629
## NA's     1840.00000  1840.00000 1840.000000 1840.000000 1840.00000000
##          layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2      layer.2
## Min.       0.0000000    0.000000    0.000000    0.000000 0.000000e+00
## 1st Qu.    0.9718253    1.666667    1.322876    1.414214 2.197058e-02
## Median     1.4142136    2.877113    2.438123    2.915476 4.184848e-02
## 3rd Qu.    2.8037673    6.009252    5.967505    7.270565 6.507355e-02
## Max.      53.2846132   85.835275   64.031242   60.055761 3.334456e-01
## NA's    1840.0000000 1840.000000 1840.000000 1840.000000 1.840000e+03

-We see that there are equal NA values in each layer We get dimension information to determine number of cells in total in each layer

sectorMovings
## class       : RasterBrick 
## dimensions  : 367, 555, 203685, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -104295.5, -102963.5, -43623.24, -42742.44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1,     layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2,     layer.2 
## min values  :  24.0000000,  27.8888889,   9.1111111,   4.5555556,  -0.5187739,   0.0000000,   0.0000000,   0.0000000,   0.0000000,   0.0000000 
## max values  : 172.1111111, 255.0000000, 189.2222222, 171.1111111,   0.7372963,  53.2846132,  85.8352751,  64.0312424,  60.0557611,   0.3334456

-We see that NA values are a very small percentage of total cells (less than one percent) of all cells#they are not visible in the plots

sectorMovings
## class       : RasterBrick 
## dimensions  : 367, 555, 203685, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -104295.5, -102963.5, -43623.24, -42742.44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1,     layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2,     layer.2 
## min values  :  24.0000000,  27.8888889,   9.1111111,   4.5555556,  -0.5187739,   0.0000000,   0.0000000,   0.0000000,   0.0000000,   0.0000000 
## max values  : 172.1111111, 255.0000000, 189.2222222, 171.1111111,   0.7372963,  53.2846132,  85.8352751,  64.0312424,  60.0557611,   0.3334456
plot(sectorMovings)

-use trim to remove NA values -from documentation: Trim (shrink) a Raster* object by removing outer rows and columns that all have the same value (e.g. NA). -Check output for NA values

sectorMovingsTrim<-trim(sectorMovings, padding=0, values=NA)

-No NA values in any layers -Plote output of trim to see if any issues

sectorMovingsTrim<-trim(sectorMovings, padding=0, values=NA)
summary(sectorMovingsTrim)
##         layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1     layer.1
## Min.       24.00000    27.88889    9.111111    4.555556 -0.51877395
## 1st Qu.    30.88889    40.44444   16.666667   10.111111 -0.21001302
## Median     34.55556    51.00000   35.555556   53.555556  0.05963125
## 3rd Qu.    40.22222    64.88889   50.888889   70.222222  0.16173571
## Max.      172.11111   255.00000  189.222222  171.111111  0.73729629
## NA's        0.00000     0.00000    0.000000    0.000000  0.00000000
##         layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2    layer.2
## Min.      0.0000000    0.000000    0.000000    0.000000 0.00000000
## 1st Qu.   0.9718253    1.666667    1.322876    1.414214 0.02197058
## Median    1.4142136    2.877113    2.438123    2.915476 0.04184848
## 3rd Qu.   2.8037673    6.009252    5.967505    7.270565 0.06507355
## Max.     53.2846132   85.835275   64.031242   60.055761 0.33344559
## NA's      0.0000000    0.000000    0.000000    0.000000 0.00000000
plot(sectorMovingsTrim)

-It looks good. The moving function may have created NA values around edges, maybe where a full three by three grid could not be made.

-determine kmeans seed values

-Normalize the ten layer stack because the range of values and variation differ between the RGB and infrared bands the NDVI. The five layers measuring the standard deviation also have a different range and variation level.

sectorMovingsTrimNorm<-raster.scale(sectorMovingsTrim)
sectorMovingsTrimNorm
## class       : RasterBrick 
## dimensions  : 365, 553, 201845, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -104293.1, -102965.9, -43620.84, -42744.84  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1, layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2, layer.2 
## min values  :           0,           0,           0,           0,       0,           0,           0,           0,           0,       0 
## max values  :           1,           1,           1,           1,       1,           1,           1,           1,           1,       1

-run kmeans on the normalized ten layer stack.

performKMeans(sectorMovingsTrim,6)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 10092250)

-Assess the clusters we have found -RAQUEL any writing about how we assess clusters (generally, or in R)

-Add the the supervised report here